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Abstract 

The MACRO detector was located in the Hall B of the Gran Sasso underground 
Laboratories under an average rock overburden of 3700 hg/cm 2 . A transition 
radiation detector composed of three identical modules, covering a total horizontal 
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area of 36 m 2 , was installed inside the empty upper part of the detector in order 
to measure the residual energy of muons. This paper presents the measurement of 
the residual energy of single and double muons crossing the apparatus. Our data 
show that double muons are more energetic than single ones. This measurement is 
performed over a standard rock depth range from 3000 to 6500hg/cm 2 . 



1 Introduction 



Underground muons are the remnants of the air showers initiated by the 
collisions of primary cosmic rays with air nuclei. These secondary cosmic 
ray muons can easily cross large amounts of matter and penetrate into 
underground laboratories. The energy spectra of underground muons depend 
on the energy spectra and composition of primary cosmic rays, on their 
interactions with air nuclei and on the muon energy loss in the rock. 

In this paper the measurement of the energy spectra of underground single and 
double muons, performed with a 36 m 2 Transition Radiation Detector (TRD), 
that was installed in the empty upper part of the MACRO detector [1], is 
presented. An analysis using part of the single muon data sample was already 
published [2]. The final sample is approximately 10 times larger and a correct 
evaluation of the systematics has been performed, thus allowing us to make 
more reliable conclusions about the spectra and the average residual energy 
of single and double muon events. The TRD sub-detector, described in sec. 3, 
uses for triggering purpose and for the measurement of the event multiplicity 
the larger area (~ 1000 m 2 ) of the streamer tube and scintillation counters 
systems. The data selection is described in sec. 4. 

To evaluate the muon energy spectra from the TRD data, we used two 
complementary approaches described in sec. 6. In the first case, an unfolding 
procedure was applied, while in the second case a parameterization for the 
underground muon energy spectrum was assumed and a best fit of the spectral 
parameters was carried out. 

To study how the energy spectra of underground muons are related to the 
primary cosmic ray spectrum and composition, a dedicated Monte Carlo 
simulation, described in sec. 7, was performed. The results of the comparisons 
between data and Monte Carlo are discussed in sec. 7.1. 
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2 Underground cosmic ray muons 



The "all-particle" flux of the primary cosmic radiation can be described by an 
inverse power law energy spectrum [3] , with differential flux given by: 

I « (i) 

where 7 w 1.7 for £ < 10 3 TeV, 7 w 2.0 for 10 3 TeV < E < 10 6 TeV and 
7 ps 1.5 for £ > 10 6 TeV. 

From the primary spectrum, it is possible to evaluate the energy spectrum at 
the Earth surface of secondary uncorrelated muons, which is given, with good 
approximation, by [3]: 

dN " ~ const ■ S~'^> ■ ( l — + °;» 54 -) (2) 



dSdfl V 1 + a£ cos 9 1 + b£ cos 9 

where £ is the muon energy at the surface, a = 1.1/115 GeV and b = 
1.1/850 GeV. The first and the second term in parenthesis of equation 2 
represent the contributions of muons from n and K decays, respectively. In 
the limit of high energies, an approximate expression of the muon surface 
energy spectrum has the simple form: 

dN 

^ = const ■ £~ a (3) 

where a = 7 + 2ps3.7. The surface muon spectral index is therefore related 
to the primary cosmic ray spectral index. 

The energy spectrum of underground muons can be derived taking into 
account the process of energy loss in the rock, which is assumed to have the 
form: 

f = -(A + W (4) 

where dh is a thin rock slab (usually in g cm -2 ), A is the contribution 
from the ionization energy loss and (5E is the contribution from radiative 
processes (bremsstrahlung, pair production and muon hadroproduction). The 
parameters A and (3 are functions of the muon energy, but for practical 
purposes can be assumed as constants [3]. The quantity e = \/ j3 is called the 
critical energy and is defined as the energy value above which the radiative 
processes become dominant. 

With the above assumptions, the general solution of equation 4 is: 

E, = (£ + e) e-? h - e. (5) 
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where is the muon energy after crossing the rock slant depth h (g cmT 2 \ 
The underground muon energy spectrum can be thus obtained from equations 
3 and 5 using the following relation: 



and is given by: 



dN 
dE,. 



dN 
~d£ 



dE 



dE 

J £=£(E„) ft 



(6) 



dN 
dE,. 



const ■ (X + e(l - e~ ph )) ° . 



(7) 



From eq. 7, the average underground muon energy at depth h is: 

W - (8) 

and its asymptotic value is e/(a — 2). At great depths h, the underground 
muon energy spectrum given by eq. 7 is almost flat for E <C (E^), and then 
decreases with energy. 



3 The MACRO TRD 



Transition Radiation (TR) is the process of the emission of X-ray photons 
occurring when an ultrarelativistic charged particle crosses the boundary 
between two materials with different dielectric constants. The most important 
features of TR are that the TR yield is roughly proportional to the Lorentz 
factor 7 of the radiating particle over a wide range of 7, and the emission 
probability of a TR X-ray is of the order of a ~ 1/137. If the rest mass m 
of the radiating particle is known, a measurement of its Lorentz factor 7 also 
allows one to evaluate the energy as E = m c 2/ y. TRDs can provide energy 
measurements over ranges typically spanning one order of magnitude. 

Due to the characteristic dependence of the TR yield on the Lorentz 
factor 7, TRDs were proposed [4,5] for the measurement of energies of 
underground cosmic ray muons in the TeV region. The TRD operated inside 
the MACRO [1] detector in the Gran Sasso underground Laboratory (LNGS) 
collected data from April 1995 to June 2000. It was designed to be sensitive 
to the energy region between 100 GeV and 1 TeV. 

The detector consisted of three modules covering a total horizontal area of 
~ 36 m 2 . Each module was composed by eleven 10 cm thick radiator layers, 
interleaved by ten planes of 32 proportional tubes 6 m long and with a square 
cross section of 6 x 6 cm 2 . A detailed description of the detector is given in [2,5]. 
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The radiator material (Ethafoam 220) contains cells of ~ 35 /im wall thickness 
and ~ 900 /im spacing, ensuring a threshold Lorentz factor 7^ pa 10 3 and a 
saturation Lorentz factor r y sat pa 10 4 , that correspond to the muon energy- 
range between ~ 100 GeV and ~ 1 TeV. 

The proportional tubes were filled with an Ar(90%) — C 02(10%) gas mixture 
and were operated at a gain of ~ 10 3 . They were equipped with a cluster 
counting read-out electronics. Wire signals were sharply differentiated and 
compared to a threshold corresponding to an energy deposit of ~ 5keV. 
In this way it was possible to discriminate <5-ray background from X-ray 
photoelectrons producing pulses exceeding the threshold amplitude. For each 
event the pulses with amplitude greater than the threshold ("hits") were 
counted in all the proportional tubes. 

The third TRD module was partially equipped with ten 1 mm thick aluminum 
foils of 2 x 2 m 2 area, that were inserted between each radiator and the tube 
plane below, in the terminal part of the module. The aluminum foils absorbed 
the TR emitted by muons in the upstream radiator layers. Using this technique 
a sample of muons were collected with only the ionization loss measured. 

A reduced scale prototype of the TRD was exposed to a pion/electron beam [5] 
to evaluate the detector response function. The physical observable related to 
the energy of the particles crossing the detector is the total number of hits 
produced in the proportional tubes. For a sample of particles crossing the 
detector with a fixed energy and at a fixed angle, the number of hits are roughly 
Poisson distributed, with an average value of few units, which depends on the 
beam energy and crossing angle. In Fig. 1 the average number of measured 
hits in the proportional tubes is plotted as a function of the Lorentz factor 7 
of the incident particles, for several beam crossing angles. Below the threshold 
value (7 = 10 3 ) there is only the contribution of the ionization energy loss; for 
10 3 < 7 < 10 4 , the TR contribution is also present and the average number 
of hits increases logarithmically with 7. 



4 Data selection and analysis 

We considered the data collected by all the three MACRO TRD modules 
during the acquisition period from April 1995 to December 2000. Two classes 
of events were analyzed: 

(1) "single muons", i.e. events with one muon in MACRO crossing one of the 
TRD modules; 

(2) "double muons", i.e. events with two muons in MACRO and only one 
muon crossing one of the TRD modules, like the one shown in figure 2. 
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To evaluate the muon energies, we associated to each muon track the hits 
produced in the TRD proportional tubes. The muon was tracked with the 
standard MACRO procedure [1], which uses the information of the streamer 
tube system. The distribution of the distances of the TRD hits from the 
expected position, calculated with the muon reconstructed track, has a 
gaussian shape with a standard deviation o ps 2 cm. We associated to each 
muon all the proportional tubes hits within 3<r from the track. To avoid badly 
reconstructed tracks, only muons crossing at least three layers of streamer 
tubes in the lower part of MACRO were used. Muons accompanied with 
electromagnetic showers initiated in the rock surrounding the detector were 
also discarded. Since the TRD was calibrated (see Fig. 1) with particles 
crossing all the ten layers of proportional tubes at angles below 45°, in the 
present analysis only muons crossing the whole detector with zenith angle 
smaller than 45° were included. 

Runs in which the TRD modules were affected by stability problems or 
malfunctioning were discarded. Also runs with reconstructed muon rates 
differing more than three standard deviations from the average values were 
disregarded, as well as those runs whose duration was less than 1 h. The final 
data samples consist of 250290 single muons and 17942 double muons, for a 
total life time of 2586 days (see Table 1). 

Fig. 3 shows the distributions of the number of hits produced in the TRD 
proportional tubes along the muon tracks for the final data samples. Since 
the second and third TRD modules were equipped with a different read-out 
electronics from that of the first module, two different response functions 
(described in the next section) were necessary to analyze the data samples 
from the first module and those from the second and third modules. The 
results of the two analysis were then combined. 



5 The TRD response function 

5. 1 Evaluation of the TRD response function 

We define N(k \ h, 6) as the distributions of the TRD hits for a sample of 
muons with zenith angle 9 that crossed a rock thickness h. The rock thickness 
h was calculated from the direction (8, <p) using the Gran Sasso map [6] and 
was converted into standard rock according the prescriptions of [7]. These 
distributions can be related to the residual muon energy spectra N(E | h, 9) 
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by: 



N(k \h,9)=^2 p(k | E, 9)N{E \ h, 9) (9) 

E 



where p(k \ E, 9) is the TR detector response function, i.e. the probability to 
observe k hits along the track of a muon with underground energy E = E^, 
crossing the detector at a zenith angle 9. 

To reconstruct the muon energy spectrum N(E \ h, 9) starting from the 
measured hit distributions (Fig. 3), once the TRD response function p(k \ E, 9) 
is known, eq. 9 must be inverted. We derived two matrices p(k \ E, 9) (one 
for the first TR module and another for the second and third modules) on the 
basis of the calibration data taken by a reduced scale prototype exposed to a 
pion-electron test beam at CERN [5]. 

We simulated, with a full GEANT-based [8] Monte Carlo, a muon sample 
distributed according a flat energy and solid angle spectrum (i.e. f^ N dn = 
constant). On the basis of the calibration data shown in Fig. 1, seven energy 
bins and four angular bins were defined. The first energy bin covers the range 
from to 50 GeV (i.e. < 7 < 500), and the last one covers the range above 
the TRD saturation (E > 1 TeV, i.e. 7 > 10 4 ). 

The number of hits produced in the TRD by a simulated muon of energy E 
crossing the detector at zenith angle 9 was extracted by the corresponding 
calibration data set. The TRD response function (a 31 x 7 x 4 matrix) was 
calculated as: 

where M{k \ E, 9) is the number of simulated muons with residual energy E 
and zenith angle 9, producing k hits, with < k < 30 and N{E, 9) is the 
total number of simulated muons with energy E crossing the TRD at a zenith 
angle 9. The simulated data had the same format of experimental data and 
were processed by the same analysis tools used for real data. 



5.2 Check of the TRD response function 



The accuracy and the time stability of the TRD response function also 
needed to be taken into account. Although the nominal gas gain of the TRD 
proportional tubes was the same of the tubes used in the prototype, during the 
long time-scale data acquisition period their operating conditions were affected 
by fluctuations. This was due to drifts of some important parameters, like the 
gas mixture composition, the atmospheric pressure and the temperature. 
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A first check of the TRD response function for energies below the threshold 
for the emission of TR was made using stopping muons, i.e. muons that are 
absorbed in the lower part of MACRO, that can be easily tagged by imposing 
simple geometrical cuts. Since the average residual energy of these muons 
is below a few GeV, we simulated a sample of muons with a flat energy 
spectrum up to 10 GeV. The same algorithm for the selection of stopping 
muons was applied to both the real and the simulated data samples. The 
measured hit distributions of stopping muons are well reproduced by the 
Monte Carlo simulation (see Table 2), thus confirming the reliability of the 
TRD calibration in the low energy region. 

A further check of the detector response function was performed using the 
data from the part of the third module equipped with aluminum foils. Since 
the aluminum foils absorbed the TR produced in the upstream radiator, this 
data sample allows to check the detector response for muons releasing energy 
in the proportional tubes only by ionization. To take into account the weak 
dependence of the TRD response on ionization (see Fig. 1), we simulated a 
sample of muons with the characteristic energy spectrum at MACRO depth 
(eq. 7) and with the same angular distribution as real data. In Fig. 4 the values 
of the average number of hits along the muon track versus the muon zenith 
angle from the Monte Carlo and from the real data are shown. The x 2 / D.o.F. 
value (evaluated using only statistical errors) is 3.3/9, showing that there is a 
good agreement between the data and the simulation. 

From this analysis we estimated the systematic uncertainty arising from 
the fluctuations of the detector operating conditions. For each of the four 
angular bins used for evaluating the TRD response function, we compared the 
measured hit distributions N(k \ 9) with the simulated ones. We remarked 
that the differences between the experimental average values and the Monte 
Carlo predictions are within ±5%. For this reason, to take into account the 
fluctuations of the detector operating conditions, we associated to the average 
values of each hit distribution of the calibration data set, an uncertainty of 
±5%. These uncertainties will propagate to the systematic errors affecting the 
measurement of the residual energy of muons. 

The systematic errors quoted in the present analysis (see sec. 6) are greater 
than the ones quoted in our previous work [2]. When that analysis was 
performed, the TRD module equipped with aluminum foils was still not taking 
data. For this reason, we could estimate the systematic uncertainties only on 
the basis of the beam test data and we decided to associate to the calibration 
data set an uncertainty of ±2%, that is lower than the value we adopted for 
this analysis. 
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6 Evaluation of the muon energy spectra 



As shown in section 5, the residual energy spectra of underground muons are 
related to the measured hit distributions in the TRD by eq. 9. 

Two different methods were applied to reconstruct the underground muon 
energy spectra by inversion of eq. 9. The first approach uses the same unfolding 
procedure described in [2] . In the second approach a parameterization for the 
muon energy spectrum was assumed and the parameters were derived using a 
best fit procedure. 

6. 1 The unfolding procedure 

Unfolding techniques are widely applied in problems where matrix inversion 
is required [9]. To reconstruct the energy spectra of single and double muons 
starting from the hit distributions of Fig. 3, we applied an iterative unfolding 
procedure [9,10]. As an initial energy spectrum (used as a starting point 
for the unfolding procedure), we assumed eq. 7 for both the single and the 
double muon events, with parameters: a = 3.7, f3 = 0.383 10~ 3 hg _1 cm 2 and 
e = 620 GeV. The final reconstructed spectrum does not depend on the choice 
of the initial spectrum. It only affects the time needed for the procedure to 
converge [9,10]. 

The energy spectrum reconstructed after each iteration is used as input for the 
next iteration. At the end of each iteration a x 2 test is performed between the 
reconstructed energy spectra and the input energy spectra. The x 2 is defined 
as: 



where iV n _i(£' i hj,9k) and N n (Ei hj,9k) are the energy spectra 
reconstructed after the (n — l)-th and the n-th iteration and a n -i(Ei | hj,9 k ), 
while cr n (Ei | hj,9 k ) are the associated errors. The iterative procedure stops 
when the energy distribution reconstructed after the n th iteration has a x 2 
probability greater than 99% to be equivalent to the one reconstructed after 
the (n— l) th iteration. The unfolding procedure was separately applied to the 
data samples from the first and from the second and third TRD modules and 
the results were combined. 

Fig. 5 shows the reconstructed differential energy spectra of single and 
double muons at the Gran Sasso underground laboratory depth. The error 
bars in the figure were calculated by adding in quadrature statistical and 



x 2 = E 



i,3,k 



[N^Ejlh^-N^jEilh^ek)] 2 



(ii) 
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systematic errors. The systematic errors are originated by fluctuations of 
the TRD response function, and were estimated by a +5% and a —5% 
variation of the calibration data. The unfolding procedure was replicated using 
the "perturbed" detector response functions; the systematic uncertainty was 
evaluated as the difference between the two results. 

The unfolding procedure reconstructs the shapes of the spectra up to 1 TeV. 
For energies greater than 1 TeV, where the TRD response is saturated, the 
spectral shapes cannot be reconstructed, and only the number of events can 
be evaluated. The average value of the energy of underground muons with 
energy above threshold, (E) E>Eo , was calculated from eq. 7 as: 



The average muon energy was evaluated as: 

(E) = (l-f)(E) E<Eo + f(E) 

E>E (13) 

where (E) E<Eo is the average energy of muons with energies up to E — 1 TeV, 
/ and (E) E>Eo are the fraction and the average energy of muons with energies 
greater than 1 TeV, respectively. For a 3% variation of the parameters a, (3 and 
e, as is typically quoted by various authors (e.g. [11]), the uncertainties on the 
average muon energies are less than 1% and are significantly smaller than our 
quoted errors. The results are shown in Table 3. The average energy of single 
muons and double muons are respectively of 270±3(stat.)±l8(syst.) GeV and 
381 ± 13(stat.) ±21(syst.) GeV. The single muon result is not in contradiction 
with the result of our previous analysis [2], where the systematic effects were 
underestimated. 



6. 2 Best fit of the spectral indices 



An alternative way for evaluating the muon energy spectra from the measured 
hit distributions is that of assuming eq. 7 as an analytic description for the 
spectra, and to derive the parameters using a best-fit. Substituting the trial 
energy spectrum into eq. 9 and summing over the zenith angles we get: 



N(k | h) = J2T,P( k I E ' d ) C (M) [E + e(l 
e e 



-/3h 



(14) 



where the normalization constants c(h, 6) represent the number of muons 
detected in each bin of depth and zenith angle. We fixed the values of the 
(3 and e parameters to f3 — 0.383 10~ 3 hg _1 cm 2 and e = 620 GeV [12]; a was 
left as the free parameter for the single and double muon data. 
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For each value of a in the range from 2 to 6 and step 0.1 we built a set of hit 
distributions N(k \ h) according eq. 14, and for each set of distributions we 
evaluated the x 2 as: 

^^ [N(h\h)-N(h\h)f 
X = 2^ 2^2^ o, u | i \ , i M (15) 



TRD modules 



a 2 (k | h) + a 2 (k | h) 



where N(k | h) is the measured hit distribution for muons (single or double) 
crossing a rock slant depth h; a(k \ h) and a(k \ h) are the Poissonian errors 
on N(k | h) and N(k \ h), respectively. 

The curves representing the x 2 as a function of a are continuous both for single 
and double muons and show a well defined minimum. The value of the best 
fit spectral index for single muons is ct\ = 3.79 ± 0.02(stat) ± O.ll(syst), with 
Xminl 'd.o.f '. = 1.51, while for the double muons is a<i = 3.25 ± 0.06(stat) ± 
0.07(syst), with x m i n / 'd.o.f '. = 0.55. The statistical error on a was evaluated 
using the values and corresponding to x 2 — Xmin + 1- The systematic 
error was estimated from the positions of the new minima of a obtained 
using the two TRD response functions evaluated from the sets of "perturbed" 
calibration data. 

The result obtained for the single muon spectral index is consistent with that 
obtained from the MACRO measurement of the underground muon intensity 
as a function of the rock depth [13]. It is also consistent with the results of 
the NUSEX experiment [14], that found a spectral index of a = 3.9llQ;3g for 
a sample of events mainly composed by single muons. 



6.3 Comparison between the results 



Fig. 6 shows a comparison between the differential muon energy spectra 
reconstructed with the two methods described in sections 6.1 and 6.2 both 
for single and double muons. In the case of single muons we find a value of 
X 2 /d.o.f . = 1.3/7, while in the case of double muons we find x 2 /d.o.f . = 8.5/7. 
In the case of single muon data, there is a good agreement between the result 
obtained using the unfolding procedure and that obtained using the best 
fit method. In the case of double muon data, there are some discrepancies, 
especially for low residual energies. 

To understand these difference, it should be noted that the result of the 
unfolding procedure does not depend on the trial spectrum. The spectrum 
used for the fit (eq. 7) is derived assuming that the surface muon energy 
spectra obey a power law with spectral index a, and in the muon propagation 
formula, e and (3 are constant, i.e. the cross sections for the radiative processes 
do not depend on the muon energy. 
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If the correct expression (eq. 2) is assumed for the surface energy spectrum 
and the dependence of e and (3 on the muon energy is taken into account, it 
is impossible to obtain an analytical expression for the underground muon 
spectrum. Hence, eq. 7 represents only a useful parameterization for the 
muon underground spectra, whose accuracy is limited by the hypotheses from 
which it has been derived. The difference between our results obtained using 
the two techniques mentioned in sections 6.1 and 6.2 are ascribed to these 
approximations. We can conclude that the parameterization is not completely 
suitable for describing the underground muon energy spectra, especially for 
events with large underground multiplicity. 



7 Monte Carlo simulation 



As shown in section 2, the underground muon energy spectrum is related 
to the primary cosmic ray spectrum. In order to investigate the relationship 
between the underground muon energy and the primary cosmic ray spectrum 
and composition, we performed a Monte Carlo simulation using different 
composition models for primary cosmic rays. 

The interactions of cosmic rays in the atmosphere were simulated with the 
HEM AS code [15], while the process of muon propagation in the rock was 
simulated using the PROPMU code [12]. Two extreme composition models 
were assumed for primary cosmic rays: 

a. the "light model" [16], i.e. a proton-rich composition model; 

b. the "heavy model" [17], i.e. a Fe-rich composition model. 

These models assume that cosmic rays are composed by five main mass 
groups (p,He,CNO,Mg and Fe). The energy spectrum of each component 
is described by means of power laws given by: 



<m 

dE 



^£-(71+1) E < Ecut 

(16) 

^ 2 £-(72+D E > Ecut 



In eq. 16 there are 5 parameters for each primary component (i.e. the 
normalization constants K x and K 2 , the spectral indices 71 and 72 and 
the cutoff energy E cut ). These parameters are not independent: usually the 
constant K 2 is expressed as a function of the others by imposing the continuity 
of the function dNi/dE at E = E cut . The all-particle spectrum is evaluated 
by adding the contributions from all the mass groups. In Table 4 the values of 
all the parameters are summarized for each component of the primary cosmic 
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rays in both the composition models we adopted. The light model is different 
from the heavy model because there is an extra-component of protons. 

Fig. 7 shows the predictions obtained using the light and heavy composition 
models for the average energy of the underground muons as a function of the 
muon multiplicity. The gap between the predictions from the two extreme 
models increases with the muon multiplicity. The two dark points represent 
our measured values, averaged on the whole rock depth range from 3000 to 
6500 hg/cm 2 . 

To explain the average energy behavior, in Fig. 8 we plotted, for the two 
composition models and as a function of the underground muon multiplicity, 
for the p, He and Fe components: a) the average energies of underground 
muons; b) the relative contribution of these three mass groups and c) the 
average energies of the parent cosmic rays to events for a fixed underground 
muon multiplicity. In a) and c) the values for the all-particle spectrum are 
also shown. From this figure the following conclusions can be drawn: 

a. single muons are originated mainly by primary protons; 

b. the contributions of primary heavy elements become relevant when the 
muon multiplicities increase, while the contribution of protons decreases; 

c. in the "heavy" composition model the contribution of Fe nuclei is relevant 
from low multiplicities (it is ~ 15% for iV M = 2 and tends to 100% for high 
multiplicities); 

d. in the light model, the main contribution is always that of protons; 

e. for a fixed muon multiplicity, muons originated by light primaries are more 
energetic than muons originated by heavy primaries. 

From the simulation we get that the average energies of underground 
muons produced by each component are almost independent on the primary 
composition model. This is due to the fact that the differences between the 
spectral indices of the primary components in the models we are examining 
are small for energies below the "knee", as can be seen from table 4. As a 
consequence, the muons produced e.g. by He nuclei have the same average 
energy, independently on the choice of the composition model. Hence, the 
differences between the predictions of the cosmic ray composition models are 
mainly due to differences in the weights of the various components in each 
model and not to differences in their spectra. 

From Fig. 8c, the primary energy needed to produce one underground muon at 
the MACRO depth is smaller than 100 TeV, while the primary energy needed 
to produce two underground muons is about 300 — 400 TeV. This means that 
an analysis of events with only one or two underground muons is sensitive to 
the energy region below the knee of the primary spectrum. 

Our measurements were compared with the predictions from the simulation. 



14 



Fig. 7 shows a comparison between the results obtained applying the the 
analysis technique described in sec. 6.1 and the Monte Carlo predictions 
concerning the average energies of single and double muons, while in Fig. 9 
the average energies of single and double muons are plotted as a function of 
the rock thickness crossed by muons. 



7. 1 Discussion of the results 

Single muon data are not in contradiction with the predictions from both the 
composition models. In fact, since the energy gaps between the predictions 
from the two extreme composition models are of the same order of magnitude 
as the error bars, single muons do not allow one to discriminate between 
different cosmic ray composition models. The sensitivity of this measurement 
to the primary cosmic ray composition is strongly limited by systematic errors 
associated to the TRD response function. 

Also in the case of double muons, our data do not allow to perform a 
discrimination between the cosmic ray composition models. In this case, 
although the errors associated to the average muon energies are smaller than 
the gap between the predictions from the two extreme composition models, 
experimental data do not provide a clean signature in favour of a given 
composition model. 

Another comparison between the TRD data and the Monte Carlo predictions 
can be done by looking at the spectral indices. For each composition model 
we fitted the underground energy spectra of single and double muons with 
the formula 7, assuming the same values of e and (3 as in sec. 6.2 and we 
compared the fit results with the TRD data. Table 5 shows the comparison 
of the data with the Monte Carlo predictions. As in the previous case, the 
single and double muon spectral indices do not allow to perform a study of 
the cosmic ray composition. 

The sensitivity of our measurement to the primary cosmic ray composition 
is mainly limited by its precision and by the poor statistics of the high 
multiplicity muon events, that does not allow to reconstruct their spectra. 
A detector with a larger area than the MACRO TRD and with a reduced 
systematics could enhance the precision of the measurement of the single 
and double muon energies and also allow to measure the energies of high 
multiplicity muons, that are more sensitive to the cosmic ray composition in 
the energy region of the knee. 



15 



8 Conclusions 



The MACRO TRD allowed the measurement of the energies of muons 
penetrating in the Gran Sasso underground laboratories, in the standard 
rock depth range from 3000 to 6500hg/cm 2 . For reconstructing the muon 
energy spectra we used an unfolding method and a best fit procedure. The 
average energies of single and double muons, evaluated with the unfolding 
technique described in sec. 6.1, are (E) 1 = 270 ± 3(stat) ± 18(syst) GeV and 
(E) 2 = 381±13(stat)±21(syst) GeV. The spectral indices for the approximate 
parameterization of the muon energy spectra (eq. 7), evaluated with the best 
fit procedure described in sec. 6.2 are a\ = 3.79 ± 0.02(stat) ± 0.11(syst) and 
a 2 = 3.25 ± O.OQ(stat) ± 0.07(syst). 

We also performed a Monte Carlo simulation to study how the muon energy 
spectra depend on the primary cosmic ray energy spectra and composition, on 
their interactions with the atmosphere and on the muon propagation in the 
rock. Our data show that double muons are more energetic than single ones 
in the standard rock depth range from 3000 to 6500 hg/cm 2 , as predicted by 
our Monte Carlo simulation. The Monte Carlo simulation also showed that 
a measurement of the energies of underground high multiplicity muons could 
provide useful information about the primary cosmic ray composition. In fact, 
we remarked that the differences between the predictions from the various 
composition models become relevant at high muon multiplicities, where the 
contribution from heavy primaries is more significant. 

The energy spectra of single and double muons reconstructed from our data are 
in agreement with the Monte Carlo predictions obtained assuming two extreme 
cosmic ray composition models. This measurement do not allow to perform 
a cosmic ray composition study because the errors are compatible with the 
gaps between the predictions from the two models. However, a measurement 
of the energy spectra of high multiplicity underground muons could provide a 
useful tool for investigating the primary cosmic ray composition in the energy 
region above the knee of the spectrum. 
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Table 1 

Summary of the data collected by the MACRO TRD. Third module data refers 
only to the part of the module not equipped with aluminum foils. 





TRD DATA 


MC SIMULATION 




Number of 


Average 


Number of 


Average 




events 


number of hits 


events 


number of hits 


Module 1 


223 


2.06 ±0.12 


3273 


2.12 ±0.03 


Modules 2 + 3 


283 


2.05 ±0.13 


3273 


2.07 ± 0.04 



'able 2 

Stopping muon data: the measured hit distributions are compared with the 
predictions from the simulation. 





(E) E <Eo 


f 


(E) 




(GeV) 


% 


(GeV) 


Single muons 


195 ± 2 sta ± 15 SJ/S 


4.5±0.1 sta ±0.7 SJ/s 


270 ± 3 s ta ± 18 SJ/S 


Double muons 


234 ± ll sta ± 18 SJ/S 


9.0 ± 0.5 sta ± lS sys 


381 ± 13 s t a ± 21 sys 



Table 3 

Average energy of underground single and double muons, with residual energy below 
the saturation threshold (Eq = 1 TeV) of our TRD (column 2); fraction of events 
with energy above the threshold (column 3); average energy of all events below and 
above threshold (column 4). 
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Light model 
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Mg 
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Table 4 

Parameters of the primary cosmic ray energy spectra according the "light" and 
"heavy" composition models. 





Single muons 


Double muons 




a 


X 2 /d.o.f. 


a 


X 2 /d.o.f. 


TRD data 


3.79 ±0.02 ±0.11 


1.51 


3.25 ± 0.06 ± 0.07 


0.55 


Light model 


3.70 


1.55 


3.14 


0.92 


Heavy model 


3.84 


1.73 


3.36 


1.02 



Table 5 



Results of the fits of underground muon energy spectra with the formula 7. TRD 
data (with the statistical and systematic error) are compared with the Monte Carlo 
predictions obtained assuming two different composition models for primary cosmic 
rays. 
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Figure 1. Average number of hits as a function of the Lorentz factor 7 for several 
beam crossing angles. The dashed lines are drawn to guide the eye. 
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Figure 3. Hit distributions N(k) vs. k for the final single and double muon data 
samples. Since the read-out electronics of the first TRD module was different from 
that of the second and third modules, the data data samples from the first module 
and the ones from the second and third modules were analyzed separately. 
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Figure 4. Check of the TRD response function by means of muons crossing the 
region of the third module equipped with aluminum foils. Data (•) are compared 
with the Monte Carlo predictions (dashed line). 
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Figure 5. Reconstructed differential energy spectra of single and double muons. 
Statistical and systematic errors have been added in quadrature. 
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Figure 6. Reconstructed differential energy spectra of single and double muons. 
Results from the unfolding procedure (black dots) are compared with those from 
the best fit procedure (open circles). Statistical and systematic errors have been 
added in quadrature. 
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Figure 8. Physical quantities related to the underground muons produced by the 
main components of the primary cosmic rays (□ = p, A= He, = Fe, • = all 
particles) plotted as a function of the muon multiplicity: a) average energies of 
underground muons; b) fraction of underground muons; c) average energies of the 
parent cosmic rays. Light (left panels) and heavy (right panels) composition models 
were separately considered. 
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Figure 9. Average energies of single and double muons as a function of the rock 
depth. Statistical and systematic errors have been added in quadrature. The 
horizontal bars represent the width of the h bins, while the central value of each 
bin corresponds to the average value of h for that bin. The last bin extends up to 
6500 hg/cm 2 . Results from the unfolding procedure are compared with the Monte 
Carlo predictions. 
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